Coulomb gap in the one-particle density of states in three-dimensional systems with 

localized electrons 



A. L. Efros 

Department of Physics and Astronomy, University of Utah, Salt Lake City, Utah 84112 

Brian Skinner and B. I. Shklovskii 

Fine Theoretical Physics Institute, University of Minnesota, Minneapolis, Minnesota 55455 

(Dated: January 20, 2013) 

The one-particle density of states (IP-DOS) in a system with localized electron states vanishes at 
the Fermi level due to the Coulomb interaction between electrons. Derivation of the Coulomb gap 
uses stability criteria of the ground state. The simplest criterion is based on the excitonic interaction 
of an electron and a hole and leads to a quadratic IP-DOS in the three-dimensional (3D) case. In 
3D, higher stability criteria, including two or more electrons, were predicted to exponentially deplete 
the IP-DOS at energies close enough to the Fermi level. In this paper we show that there is a range 
of intermediate energies where this depletion is strongly compensated by the excitonic interaction 
between single-particle excitations, so that the crossover from quadratic to exponential behavior of 
the IP-DOS is retarded. This is one of the reasons why such exponential depletion was never seen 
in computer simulations. 



I. INTRODUCTION 

The study of localized, interacting electrons in a dis- 
ordered system has been a source of interesting physics 
for nearly half a century. The canonical example of such 
a system is a lightly doped, compensated, n-type semi- 
conductor, where electrons become localized on donor 
sites The study of electron-electron interactions in 
the localized regime was initiated by Pollak 0] and Srini- 
vasan [sj . Efros and Shklovskii [4, 5] later argued that the 
single particle density of states (IP-DOS) tends to zero 
at the Fermi level as a result of the long-range part of the 
Coulomb interaction between electrons. They proposed 
the following universal soft gap in the IP-DOS near the 
Fermi level at temperature T = 0, which is called the 
Coulomb gap and depends only on the dimensionality D 
of the system: 



g(e) = —\e\forD = 2 



and 



Ke) = -^e^ for Z? = 3. 



(1) 



(2) 



Here the reference point for the energy e is the Fermi 
level and denotes the square of the electron charge 
divided by the dielectric constant. Eqs. ([1]) and ([2]) were 
derived for the case when the bare DOS, which is the 
DOS without Coulomb interactions, has a non-zero value 
at e = 0. 

The primary method for quantitative, theoretical 
study of the Coulomb gap is through computer simula- 
tions 043) which mostly confirm the above results for 
the IP-DOS. The most important experimental manifes- 
tation of the Coulomb gap is its effect on the variable 
range hopping conductivity a, which as a consequence of 
the Coulomb gap obeys the law [ij. 



a = (Toexp[-(To/T)i/2] 



(3) 



observed in dozens of experimental works. Here ao is a 
constant and To = /Jne^/^, where ^ is the localization 
length and is a numerical constant that depends on 
the dimensionality: w 2.8 fl| while /32 ~ 6 [11 [111 . 

Eqs. (HI) and ^ are derived from the requirement that 
the ground state be stable with respect to the transfer of 
a single electron from an occupied to an empty state. In 
this introduction we will just mention that this depletion 
of the bare IP-DOS is a result of the excitonic attraction 
of electrons and holes near the Fermi level, which elimi- 
nates electrons and holes with small |e| that are close in 
space, pushing them to higher energies. In other words, 
one can say that electrons and holes that remain in the 
ground state and contribute to the IP-DOS "intimidate" 
each other. 

While Eqs. ^ and ^ reflect the effect of this first- 
order stability criterion, a more challenging question is 
how the IP-DOS is affected by higher-order stability cri- 
teria, in which two or more electrons are transferred si- 
multaneously. It has been shown that such criteria do 
not significantly modify the linear gap in the 2D case [Eq. 
(H))], but they should modify the quadratic gap in the 3D 
case [Eq. ([2])]. Previous authors [1, [13 have suggested 
that compact dipole excitations intimidate single parti- 
cle excitations, thereby making the IP-DOS exponen- 
tially small in the limit of small enough energies. It was 
also understood that this effect does not modify the law 
of variable-range hopping conductivity [Eq. ([3|)] because 
when dipole excitations play an important role hopping 
conductivity occurs by multi-electron excitations called 
electronic polarons, for which the DOS always obeys Eq. 
© [il,i,[ll- The depleted IP-DOS could in principle be 
seen in tunneling experiments [l9| . 

However, such strong, exponential depletion of the IP- 
DOS was never seen in computer experiments (^-[T^]. For 
large simulated samples this is the case because computer 
simulations are not able to reach the ground state of the 
system or cannot enforce a sufficiently large number of 
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higher-order stabihty criteria. In this paper we show that 
even for samples where the ground state can be reached 
the crossover from a quadratic to an exponentially small 
IP-DOS with decreasing energy |e| is strongly retarded. 
The reason is that the survival probability P{e) of sin- 
gle particle excitations in the presence of intimidating 
dipoles remains a weaker function of e than Eq. ^ over 
a wide range of energies. The crucial idea of this paper is 
that the bare DOS multiplied by P(e) can be considered 
as an effective new bare DOS. We show that because of 
the remarkable stability of ^(e) with respect to changes 
in the bare DOS, such modification produces only a small 
effect on the IP-DOS over a substantial range of inter- 
mediate energies. 

Finally, at even smaller energies, where the exponen- 
tial effect of dipoles reduces the modified bare DOS so 
strongly that the mutual electron-hole intimidation plays 
a weak role, the IP-DOS decreases exponentially and is 
proportional to P(e). At present, computer simulations 
are not able to find the ground state in systems that are 
large enough to explore such small energies. 

The remainder of this paper is organized as follows. 
In Sec.|TT]we present the simple model Hamiltonian that 
is used in most numerical studies of the Coulomb gap 
and we discuss the role of mutual intimidation of sin- 
gle particle excitations, including its description by the 
self-consistent equation (SCE) [5[[ll]. In Sec. |lll]we re- 
capitulate old theoretical results related to dipole excita- 
tions and their mutual intimidation We then discuss 
the depletion of the bare DOS by surviving dipole excita- 
tions and we use this depleted DOS as a bare DOS 
to solve the SCE for g{e). In Sec. |IV] we study analyt- 
ically the stability of the Coulomb gap for hypothetical 
cases where the bare DOS follows a power law, and we 
confirm all qualitative conclusions obtained numerically 
in Sec. HIl 



II. EXCITONIC EFFECT ON SINGLE 
PARTICLE DENSITY OF STATES 

In order to remind the reader of the derivation of Eqs. 
([!]) and we start from the model Hamiltonian 



H 



E'^^"» + ("»-l/2)K--l/2). (4) 



The electrons described by this Hamiltonian are assumed 
to occupy sites on a cubic lattice with lattice constant /. 
ni = 0, 1 is the occupation number of site i and ry is the 
distance between sites i and j. The quenched random site 
energies (f>i a-re distributed uniformly within the interval 
[—Ae^ /I, Ae^ /I], where ^ ^ 1 is some positive number, 
so that the bare DOS created by the random site energies 
is the "box function" goi^) — goo0[^ — kl/(e^/0]- Here, 
goo — {2Ae'^P)~^ and Q[x] denotes the Heaviside step 
function. In order to make the system neutral each site is 
given a positive background charge e/2. Due to electron- 
hole symmetry, the chemical potential ^ = 0. The single 



particle energy at site i is given by 

=2 



1/2). 



(5) 



At zero temperature (in the ground state of the sys- 
tem) all states with < are occupied and all states 
with > are empty, since the ground state must be 
stable with respect to transfer of an electron from an oc- 
cupied site to infinity or from infinity to an empty site. 
This condition is only the first stability criterion of the 
ground state. The second stability criterion considers the 
transfer of one electron from some site i that is occupied 
in the ground state to a site j that is vacant in the ground 
state. The change in H resulting from such a transfer. 



^ij ~ „ ' 



(6) 



must be positive for the stability of the ground state. 
The last term in Eq. ^ reflects the simple fact that 
the ground state energy ej includes the potential of the 
transferred electron, which initially was at site i. 

One can see the origin of the Coulomb gap from Eq. 
Since ej > and < 0, the first two terms give a 
positive contribution to f2y , while the third term is neg- 
ative. Thus, any two sites i,j with energy close to the 
Fermi level should be separated in space to keep ilij > 0. 
To see how this produces the Coulomb gap, consider sites 
whose energies fall in a narrow band (—e/2, e/2) around 
the Fermi level. According to Eq. any two sites in 
this band with energies on opposite sides of the Fermi 
level must be separated by a distance not less than 
e^/e. Therefore the concentration of such sites n{e) can- 
not exceed (e/e^)^. Thus, the IP-DOS 5(e) = dn{e)/de 
must vanish when e tends to zero at least as fast as e^~^. 

In this way we arrive at a simple upper bound for the 
IP-DOS. Applying the principle of the micro-canonical 
distribution to the disordered system we find that all 
states except the forbidden ones are homogeneously dis- 
tributed in phase space. Thus, the bound g{e) cx e^~^ 
is an exact one if we enforce only the criterion Qij > 0. 
Applying all other stability criteria can only reduce the 
IP-DOS. 

The numerical coefficients in Eqs. ([T]) and ([2]) were 
later found by solving the SCE for g{e), which is de- 
signed to describe mutual intimidation of single particle 
excitations. In 3D this equation has the form [sl. Il8j| 



5(e) = 50(e) exp 



2tt 



(|e|+e')3 



(7) 



where 50(e) is the bare DOS. 

The asymptotic solution of Eq. Q at small energies 
e < A leads to Eq. ^ (5|]. Here A is the width of the 
Coulomb gap, defined as the energy e at which Eq. ([2]) 
matches 500, so that A = {e^ /l){Tr/6Ay^^ < A{e'^/l). 
Solution of the 2D version of the SCE at small energies 
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e <C A leads to Eq. ([T]). Remarkably, these asymptotic 
solutions do not depend on A or Z and in this sense are 
universal. It should be noted that Eq. ([7|) does not de- 
scribe the corresponding increase in high energy states 
that accompanies the depletion of low energy states. In 
this sense its solution g(e) is not exactly normalized to 
contain the same total number of states as goi^e), but this 
should not be very important for A 3> 1 . 

At only moderately small energies, where e/A is not 
asymptotically small, the integral equation of Eq. (O has 
no known analytical solution. Nonetheless, one can find 
numerical solutions to the SCE by using a simple itera- 
tion procedure, wherein an initial guess is made for the 
function g{e) [say, that of Eq. ([2|)] and then inserted into 
the right-hand side of Eq. (O to produce a new estimate 
for (?(e). This iteration can be repeated an arbitrary num- 
ber of times until a convergent solution is found for the 
hmction g{e) over the desired range of energies [2d| . 

The result of this procedure is shown for ^ = 2 in 
Fig. [T] by the thin red line. At small energies this solu- 
tion approaches its asymptotic form [Eq. shown by 
the dotted black line] . In the dimensionless units of this 
plot the solution ioi A = A (not shown) coincides almost 
exactly with the solution for ^ = 2. 



III. DIPOLE EXCITATIONS AND THEIR 
INTERACTION WITH SINGLE PARTICLES 
EXCITATIONS 

Until now we have taken into account only the condi- 
tions that H be minimized with respect to the change of 
one or two occupation numbers. An analysis of the role 
of stability criteria with respect to the change of three 
or more occupation numbers shows that in the 2D case 
the single-electron approach, used above, is good, while 
in the 3D case the physics is more complicated 

It is convenient to talk about this problem by intro- 
ducing a dipole branch of excitations. For a given elec- 
tron transition from a site that is occupied in the ground 
state to an empty site, the energy is given by Eq. 1^. 
The result of such a transition is the formation of an 
"electron-hole pair". If the energy $7^ of this transi- 
tion satisfies flij » e'^/rij, then the pair can be thought 
of as two independent single-particle excitations: an ex- 
tra electron on site j and an extra hole on site i. If 
^ij <C e^/rij, on the other hand, the pair constitutes a 
strongly bound small dipole [1, [3, [i3| ■ The typical arm 
of such dipoles is tq = e^/A. At small O these dipoles 
are spatially separated from each other and do not par- 
ticipate in dc conductivity. They do contribute, however, 
to the ac conductivity and the low temperature thermo- 
dynamics [i. [sl. [iTj . 

In the first approximation the interaction between 
these dipole excitations was assumed to be weak f^J^ and 
their density of states (2P-D0S) was predicted to 

be constant at small energies and given by $(^2) = 500- 
A more careful consideration [3 shows that dipole- 



dipole interaction leads to some mutual intimidation, re- 
sulting in the logarithmic depletion of ^'(S^). Specifically, 
<i>(r2) — goo/[lii(A/r2)]^/^. At the same time, the remain- 
ing soft dipoles have dipole arms that are shorter than ro, 
with an average arm ro/[ln(A/r2)]^/'*. This happens be- 
cause shorter dipoles interact more weakly and therefore 
survive with larger probability. 




FIG. 1. (Color online) The IP-DOS g{e) (thick black curve) 
and the probability P{e) for a single particle excitation to sur- 
vive interaction with dipole excitations (blue, dashed curve). 
The thin red curve, 51(e), is calculated from Eq. the green 
dashed-dotted curve is the product of this solution with P(e), 
and the black dotted line is the asymptotic result of Eq. ([2|. 
Both (/(e) and (71(e) are calculated using A = 2. 

Let us now describe the interaction between dipole ex- 
citations and single particle excitations. We first note 
that dipole excitations are more abundant than single 
particle excitations due to the strong depletion of the 
IP-DOS at the Fermi level. Thus, it is natural to think 
that after finding the 2P-D0S one can use it to study 
the role of dipoles in intimidating single-particle excita- 
tions. Following this logic Efros [5] noticed that in 3D a 
dipole with the proper orientation and proximity can in- 
timidate a single particle excitation. Using the constant 
2P-D0S $(r2) = gooj Efros showed that the probability 
that a single particle will survive this interaction is given 
by P(e) = exp[— (A/e)^/^]. An improved understanding 
of the effects of dipole-dipole interactions [l^ leads to 
a logarithmic increase of this probability. A calculation 
along the lines of Ref. [13 shows that a more accurate 
expression for the survival probability is 



Pi(e) = exp 



A 



1/2 



1/2 in-V8 - 



A 



(8) 



where 7 = 1.5 is a numerical constant JJ,]. Both the log- 
arithmic reduction of $(^2) and the logarithmic reduction 
of the size of dipoles discussed above contribute to the 
logarithmic factor in Eq. ([5]) that increases the probabil- 
ity i-i(e). The subscript 1 in Pi(e) emphasizes that it 
describes the effect of intimidation by a single dipole. 
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Later in Ref. the authors realized that coUective 
intimidation by many dipoles surrounding a one parti- 
cle excitation is stronger than intimidation by a single 
dipole. This can be seen by considering that around an 
empty site (a positive single particle excitation) there 
is a certain probability to find a cloud of many dipole 
moments oriented away from the hole. Such dipoles can 
re-polarize (move their electron away from hole) when an 
electron is brought to the empty site from a low-energy 
occupied site. This process lowers the system's total en- 
ergy. To guarantee that the original state is not vulner- 
able to such an event, and, therefore, is the true ground 
state, one has to exclude such a polarization atmosphere. 

The probability of a hole or electron surviving in the 
ground state was calculated in Ref. 



P(e) = exp 



-7 



(9) 



Remarkably, the argument of the exponential in P{e) 
is exactly the square of the argument of the exponential 
of -Pi(e). Thus, in the interesting case where the latter 
argument is larger than unity, the intimidation of single 
particle excitations is determined by the collective effect 
result of Eq. ©. 

Eq. © is of course derived in the limit that ln(A/e) ^ 
1. Attempting to use it at A/e ~ 1 leads to nonphysical 
divergence, since ln(A/e) tends to zero. Instead, the cor- 
rect behavior for P(e) is that it should approach unity at 
e/A > 1. One can model such behavior while still main- 
taining the correct small energy asymptotic form in Eq. 
(|9|) by introducing a multiplicative "crossover function" 
/(e/A) into the exponent of Eq. ([5]): 



P(e) = exp 



A 



In 



-7/4 



/( 



(10) 



The function f{x) should have the properties f{x — > 0) = 
1 and f{x > 1) = 0. One should also ensure that /(x) ap- 
proaches zero at x = 1 faster than (x — 1)^/^, so that the 
divergent behavior of Eq. ([9|) at e/A = 1 is eliminated. 
One obvious choice for /(x) is 



/(x) = (i-x)"e(i-x) 



(11) 



with r] > 7/4. In the numeric results presented in Fig. [T] 
we use 77 = 4, which produces a smooth crossover from 
the asymptotic behavior of Eq. ^ to P(e) = 1 at e > A. 

We can now consider how the IP-DOS is affected by 
the elimination of single particle states, as given by the 
small survival probability in Eq. (jlOp . Asymptotically at 
very small e, P{e) goes to zero faster than e^ and should 
play a dominating role in the IP-DOS. This is why the 
fact that numerical data for g{e) in 3D are close to Eq. 
^ looked puzzling 

We would like to argue here that the dominance of Eq. 
([T0| requires very small e. One can notice that around 
e/A = 0.3 the probability P(e) still has the strength of 
the first power of e/A. It reaches the strength of (e/A)^ 



only at e/A « 0.1. This is seen in Fig. [1] where the 
probability P(e) is plotted by the dashed blue line. 

Still, one may think that the effect of mutual intimi- 
dation of single particle excitations and the effect of in- 
timidation of single particles by dipole excitations are 
multiplicative, so that to obtain g{e) one should multiply 
the solution of Eq. ^ for the bare "box" DOS, which 
we dub gi{e) (thin red hue), by P(e) from Eq. ([TOl) . This 
product is shown by the dashed-dotted green curve in 
Fig. [II 

We argue that instead of directly multiplying the final 
result of Eq. ^ by P(e), one should think that the role of 
P(e) is to modify the bare DOS go{^)- Thus, to find g{e) 
one should replace goi^) in Eq. (O by go{€)P{e) and solve 
this new self-consistent equation for g{e). The result of 
this calculation is shown by the thick black line in Fig. 
[TJ as determined by the numerical iteration procedure 
described in Sec. |lTl 

One can see that at e/A > 0.1 the IP-DOS stays much 
closer to Eq. ([7]) than the product P(e)(;i(e). Thus, in 
the range of energies e > 0.1 A, the solution of Eq. ([7]) 
is very stable with respect to modification of the bare 
density of states. The mechanism of this stability is as 
follows. When go{e)P{e) becomes moderately small, the 
exponential factor on the right-hand side of Eq. ([T]) grows 
sharply, supporting an almost unchanged value for g{e). 
In other words, any attempt at moderate reduction of 
g{e) by dipole intimidation is met by resistance from the 
compensating effect of weakening of the intimidation by 
single-particle excitations. This is an extension of the 
universality of the small energy solution of Eq. ([7]) with 
respect to varying the bare "box" density of states goo 0, 

Returning to Fig. [U we note that at e ^ 0.1, where 
the probability P(e) changes substantially faster than 
(e/A)^, our result for 5(e) follows the exponentially small 
P(e). At present, however, this is a purely theoretical 
range, since computer simulations can find the ground 
state at A ~ 3 only for lattices of size 10 x 10 x 10 or 
smaller For such small samples, finite size effects Q 
limit the range over which one can determine the IP- 
DOS of infinite systems to e/A > 0.3 [2l|. At these 
energies the value of g{e) in the ground state follows Eq. 
@ l2l|, in agreement with our prediction in Fig. [1] (the 
thick black curve). 



IV. STABILITY OF THE COULOMB GAP: 
ANALYTICAL SOLUTIONS 

In the previous section we dealt with the stability of 
the IP-DOS Coulomb gap with respect to moderate de- 
pletion of the bare DOS go(e) by solving the SCE nu- 
merically. In this section we would like to illustrate this 
stability for an analytically solvable model. We consider 
this model only in its 3D version, but the two dimensional 
version is similar. 

Assume that the bare "box" DOS (70(e) — gooQ[A — 
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|e|/ (e^/Z)] used above is modified at small energies e < A 
by the power law factor 



50 = .goo(e/A)^ 



(12) 



where a > 0. In order to find out how this changes 
the low energy solution of Eq. ([7]), we write it in the 
differential form 



ding dingo 



de 



de 



(13) 



At small values of e and a < 2 the solution of this 
equation is 



(14) 



where /3 is a numerical factor. Substituting Eq. into 
Eq. ([Ta]) one gets 
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3(2 - a) 
2^ ' 



(15) 



At a = we are back to the bare "box" go(e) and Eq. 
Pl)) coincides with Eq. One can see that even at 

< a < 2 the density of states is still quadratic and the 
only effect of the depletion of the bare DOS in Eq. 
is to reduce the numerical coefficient (3. In this way at 
small energies the solution g{e) is remarkably stable to 
the strong change in ^q. The reason for this stability is 
that at < a < 2 the number of available sites in a small 
energy band near e = exceeds the number of electrons 
that is permitted by their interaction. Of course, the 
electrons cannot arrange themselves as comfortably as in 
the case of constant go, and this fact is reflected in the 
reduction of the coefficient (3. 

On the other hand, when a > 2 the number of places 
for electrons in a small energy band near e = is less than 



the number of electrons permitted by their interaction. 
Therefore, g{e) has the same energy dependence as go{f^)- 
This can be seen formally from Eq. Substitution of 
.go(e) into the integral term with a > 2 allows one to set 
e = inside the integral at e/A ^ 1. Then at small e 
one gets g{e) = Bgo{e), where _B is a constant given by 



B = exp 



2tt 



9o{e')de' 



(16) 



Here we have used the fact that 50(e) is zero at e > 
A{e'^/l), so that the integral in Eq. ([T6|) converges. 

In the critical case, a = 2, one can show that the 
asymptotic solution at e — > is 



.9(e) 



27re6 ln(A/e) + /i(A/e)' 



(17) 



where h{A/e) is a small correction to ln(A/e). We were 
not able to find this correction analytically, but numerical 
solutions to the SCE suggest that h{A/e) ~ 7.3 at e/A < 
10-3. 

The analytical solutions of this section confirm our nu- 
merical experience from the previous section. Namely, 
that when the bare density of states go{e) at small ener- 
gies is larger than the one given by Eq. the solution 
of Eq. ([7]) is very close to Eq. On the other hand, 
when (70(e) is smaller than Eq. ([2]) at small energies the 
density of states g{e) is close to (70(e). 
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